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ABSTRACT 

Massive stars are known to have a high multiplicity, with examples of higher order 
multiples among the nearest and best studied objects. In this paper we study hierar- 
chical multiple systems (an inner binary as a component of a wider binary) of massive 
stars in a clustered environment, in which a system with a size of 100-1000 au will 
undergo many close encounters during the short lifetime of a massive star. Using two 
types of Wbody experiment we determine the post-formation collision probabilities 
of these massive hierarchies. We find that, depending on the specifics of the environ- 
ment, the hierarchy, and the amount of time that is allowed to pass, tens of percent 
of hierarchies will experience a collision, typically between the two stars of the inner 
binary. In addition to collisions, clusters hosting a hierarchical massive system produce 
high velocity runaways at an enhanced rate. The primordial multiplicity specifics of 
massive stars appear to play a key role in the generation of these relatively small num- 
ber events in cluster simulations, complicating their use as diagnostics of a cluster's 
history. 

Key words: binaries: general - methods: n-body simulations - open clusters and 
associations: general - stellar dynamics 



1 INTRODUCTION 

Collisions between stars is a topic that, depending on its 
context, can carry with it a whiff of controversy. Several au- 
thors have studied it as a mechanism for general massive star 
formation (e.g. |Bonnell et al.|1998||Bally fc Zinnecker|2005| |, 
although the prime motivation for considering collisional for- 
mation has weakened as numerical models of star formation 
have become more sophisticated. More recent work on col- 
lisions in the formation phase of massive stars has tended 



to focus on the most massive, exotic objects (Davis et al 
|2010||Moeckel fc Clarke|2011||Baumgardt fc Klessen|2011 1 
This scenario bears many similarities to the runway colli- 
sion product route to intermediate mass black holes in core 



collapsed massive clusters (e.g. 


Portegies Zwart & McMillan 


2002 


IGiirkan et al.|2004||Preitag et al.|2006||Gaburov et al. 


2010 





lisions of stars themselves in the purely gravitational N- 
body studies devoted to it is not one of them. Given suf- 
ficiently high stellar densities and enough time, collisions 
will occur. As another example, collisions have been invoked 
in order to explain blue stragglers (e.g . |Leonard||1989 Sig 



|urdsson "fc Phinncy 1993. B acon et al.||1996[ |Fregeau et al 
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|2004[ ). While most of these studies are purely gravitational 
and involve only binaries and singles, recent work has pushed 
into higher order multiples (Leigh & Geller 20121. Beyond 



gravity, other work examines the hydrodynamics of stellar 
collisions (e 



Davies et al. 



1993 Sills et al. 



1997 



2002 



Lombardi et aL||2003| |Laycock fc S ills 2005 ; [Gaburov et al.| 



2008]|2010l. 



Our focus here is on collisions involving massive stellar 
systems (with masses greater than 10 Mq) after their for- 
mation and during their brief main sequence lifetime, when 
they are gravitationally interacting with other stars in nor- 
malilenvironments. Our main motivation is the observation 



that massive stars are highly multiple ( Garcia & Mermilliod 
|Sana et al.|[2008| |2M9l |Chmi et al.||2012| |Sana et aT 



2001 



20121; we lay the scenario we have in mind in more detail 



1.1 Encounter timescales 

Consider first encounters between a binary system and a 
passing star with a maximum closest approach between the 
intruder and the binary centre of mass r enc . For something 



1 By which wc mean non-extreme, e.g. without invoking tem- 
porarily high stellar densities at some stage of a star cluster's 
early existence. 
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Figure 1. The encounter timescale versus encounter radius for 
representative massive stellar systems. The red line shows an en- 
counter between 20 Mq of stars, the black lines for 150 Mq. Thin 
lines are in an environment with stellar number density n = 10 4 
pc -3 and one dimensional velocity dispersion u = 1 km s — 1 , and 
the thick lines are for n = 10 5 pc - 3 and <r = 10 km s — 1 . 



interesting to happen, r enc should be something like the bi- 
nary semi-major axis a. We can estimate the frequency of 
these encounters as t enc w (new) -1 , with the number den- 
sity of single stars n moving at velocity v relative to the 
binary with an encounter cross section a. The cross section 
must take into account the gravitational focussing of orbits 
(e.g. Leonardj 1989 1, so that for a binary of mass Mb and an 
intruder of mass Mi the cross section is 



1 + 



2G(M b + Mi) 



(1) 



In figure [T] we show this encounter timescal^Jas a func- 
tion of r enc for representative values of the system mass 
and cluster environment. At small encounter separations the 
gravitational focussing limit depends on the stellar masses; 
the large separation geometrical limit is set by the cluster 
environment. Even very massive systems (encounter part- 
ners totalling over 100 Mq) have encounter timescales of 
order 1 Myr for encounter radii of 10 au. At 100s to 1000s of 
au, any massive system in a typical cluster has an encounter 
timescale of a fraction of a Myr. Binaries in this separation 
range should have close encounters during the stars' main 
sequence lifetimes. 



1.2 The insufficiency of mere binaries 

Encounter frequency is not enough to cause collisions, how- 
ever; the multiple systems should also be resistant to quick 
and outright destruction by the intruder and, furthermore, 
likely to have individual stars pass closely enough together 
during the interaction to collide. A binary with component 
masses mo and mi and semi-major axis a that is energeti- 
cally 'soft' relative to the energy of a typical intruder star 



2 We actually plot the timescale averaged over a Maxwellian ve- 
locity distribution | |Binney fc Tremaine|[2008[ l, which introduces 
a correction of order unity to the nav estimate. 



of mass Mi with velocity dispersion er, i.e. 
Gmorni 



2a 



(2) 



will be quickly disrupted without undergoing any complex 
small- N dynamics ( |Heggie|[T975} |Hut fc BahcaiT|[T983l ). A 
'hard' binary by contrast, with binding energy well in ex- 
cess of the typical kinetic energy of it's neighbours, is more 
likely to undergo chaotic and complex resonant encounters. 
Massive binaries are energetically hard out to large enough 
separations to allow interesting encounters to take place on 
timescales small compared to the short lifetimes of O stars. 
However, with the wider separation relative to the stel- 



lar radius comes a decreasing chance of collisions (Fregeau 
|et al.|2004| ), and the increased encounter rate does not gen- 
erate a higher net collision frequency. Generally, for an en- 
counter to have a high probability of collisions between stars 
you need binaries that are so compact (R*/a < 10~ 3 , or 
a < 10 au for massive stars) that in a cluster with peak 
number densities of order 10 4 pc -3 the encounter rate is too 
small to be significant over a few Myr. 

What if one of the stars in a binary that is large enough 
to have frequent encounters is not a single, but itself a com- 
pact binary? Generally referred to as a hierarchical system, 
this setup introduces the possibility of a two-stage collision 
process: the wide (but still energetically hard) binary acts 
as a net to draw passing stars into encounters at a high 
rate. The interaction with the wide binary shepherds in- 
truders into close proximity with the compact binary with a 
higher frequency than it would experience by itself, poten- 



tially leading to collisions. Leigh & Sills (20111 analytically 



studied this effective enhancement of the inner binary inter- 
action rate in the context of blue straggler formation; |Geller| 
et al. (20131 numerically confirmed the predictions of that 



work. The process is relevant to dynamically formed triples 
given the time and length scales of old open and globular 
clusters, and here we investigate it in young clusters with 
massive primordial triples. 

1.3 Observational and theoretical support for 
higher order massive multiples 

Higher order multiples are found among nearby stars of 
approximately Solar mass (e.g. Duquennoy fc Mayor||199r 



Raghavan et al.||2010 |, as well as the nearby massive stars 
in the Trapezium (© Ori A and B are a triple and at least 
quadruple; Preibisch et al. 1999||Schertl et aT 2003 1 and else- 
where in Orion; O Ori A ( Preibisch et al.|2001 I and a Ori 
AB ( |Sanz-Forcada et al.|2004| . In the embedded phase there 
are already hints of high order multiplicity among massive 
objects (e.g. W3 IRS 5; Megeath et al.fl2005 ), suggesting a 
primordial rather than than dynamic origin. 

On the theoretical side, perhaps the most likely forma- 
tion scenario for primordial massive binaries is disc fragmen- 
tation (|Adams et al.||1989| |Laughlin fc Bodenheimer||1994| 



Bonnell 1994 Bonnell & Bate 1994a b). In simulations of 
star formation, starting from a variety of initial conditions 
and utilising different techniques, small-TV multiples born 



out of disc fragmentation are common (e.g. 


Krumholz et al. 


2007 


2009;|Bate 


2009 


Peters et al.|2010| |Stamatellos et al. 


2011 


Greif et al. 


2011 


Bate|2012i. Around massive stars in 



particular fragmentation seems prevalent, perhaps because 
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in contrast to Solar type stars continued high accretion rates 
onto the disc result in conditions more favourable to frag- 



bedded 7th order estimator for error control ( Prince & Dor- 



mentation into binaries of multiple systems (Kratter et al 



|2008| |2010[ > . Fragmentation in such massive discs tends to 
occur at the disc edge, at radii of order 10 2 to 10 3 au, which 
is the right separation for frequent encounters in a cluster. 

A possible scenario for the formation of a hierarchical 
triple is fragmentation followed by accretion-induced orbital 



shrinkage ( Bonnell fc Bate|2005[ ) or disc migration (Goldre 



|ich fc Tremaine||1980[ ) leading to a binary with separation 



of order 10 au, surrounded by a circumbinary disc (Arty- 
mowicz fc Lubow|1994 1996 1. Further fragmentation of this 



disc could lead to another companion at larger radii. Alter- 
natively, early protostellar dynamics involving a disc that 
fragments into multiple objects may settle into a hierar- 
chical arrangement (e.g. Sterzik & Durisen 19981, partic- 



ularly when dissipative tidal forces are included (Mardling 
fc Aarseth|2001| ) 



While hydrodynamic simulations of star formation cre- 
ate higher order multiples with some frequency, they are 
seldom included in initial conditions when performing grav- 
itational TV-body studies of star clusters (although they cer- 
tainly can form during the course of a simulation). Encoun- 
ters involving triples have only recently seen some of the 
attention given to binary-single and binary-binary encoun- 



ters (Leigh & Geller 20121, and their focus was on stars 



of roughly Solar mass. If massive hierarchies are included 
in the initial conditions, the timescale and energetics argu- 
ments above are at least suggestive of resultant interesting 
encounters. 



1.4 The plan 

In this paper we set out to numerically determine stellar col- 
lision rates involving massive, primordial, hierarchical triples 
in a clustered environment. We make use of two types of 
numerical experiment. In section [2] we describe idealised en- 
counters between isolated hierarchies and perturbing stars, 
and in section [3] we describe those results, and in section [1] 
we briefly discuss their restriction to coplanar hierarchies. In 
sections [5] and [6] we describe more expensive and thus lim- 
ited experiments involving a full cluster simulation to verify 
the small- N work. Given the history of massive star forma- 
tion theory, we would like to emphasise that the processes 
we are studying here have very little to do directly with 
the formation of individual massive stars; we are studying 
the relatively uncontroversial dynamics of star clusters us- 
ing initial conditions that can arguably be supported by the 
gamut of massive star formation theories. Our primary as- 
sumptions are the existence of star clusters and higher order 
multiples. 



2 SMALL-N EXPERIMENTAL SETUP 

For this part of the paper we set up isolated hierarchical 
triple systems of massive stars and bombarded them with 
other stars (called 'intruders'). We used the code Fewbody 
(Fregeau et al. 20041, an iV-body integrator optimised for 



this type of small- iV experiment. While Fregeau et al. (2004 1 



should be consulted for details, briefly: the code uses a vari- 
able timestep 8th order Runge-Kutta method with an em- 



mand||1981 1. During integration if a subsystem of the stars 



becomes isolated (e.g. a temporary binary is on a looping 
orbit far from the other stars) it is advanced analytically 
until the tidal force on the subsystem from other stars ex- 
ceeds some fraction of the internal force of the subsystem; 
we use 10 -5 as the tolerance for direct integration. 

The hierarchical triple system consisted of an inner bi- 
nary with semi-major axis ao which was in turn one member 
of a wider binary with ai 3> ao- Both binaries initially had 
zero eccentricity, with orbital planes randomly oriented rel- 
ative to each other. We chose the mass of the primary star 
mo from a Salpeter-like mass function dN/dm oc m -2,3 be- 
tween 10 and 150 Mq, and the inner binary partner's mass 
mi by randomly selecting a binary mass ratio q = mi /mo 
in the range 0.1 to 1.0, with the mass ratio distributed 
as a power law with p(q) oc q~ 01 . Recent observations of 



O-star binaries in six Galactic open clusters (Sana et al 



2012) motivated this choice. The outer binary partner's mass 



m,2 was likewise drawn by picking a mass ratio relative to 
the system's primary star mo- The inner binary thus had 
mass Mb = mo + mi, and the hierarchy as a whole had 
Alt = Mb + m,2. The intruding star's mass Mi had the same 
mass function as the primary, but in the range 1 to 150 Mq. 
Table [l] collects these symbols for reference. 

During integrations, Fewbody detects collisions be- 
tween stars if their surfaces are detected to have touched. 
These collisions are mass and momentum conserving, an ap- 
proach sometimes referred to as 'sticky spheres'. We therefor 
needed a radius for each star. We used a simple parameteri- 
sation of main-sequence radii: with radii and masses here in 



Solar units, r = m for m < 2.5, and r 
for m ^ 2.5 (see e.g. figure 4 in 



2.5»- t) (m/2.5) 1J 



Demircan & Kahraman 



19911. We make no claim that this choice was overly real- 



istic. Since there is no stellar evolution included in Few- 
body, we merely wanted plausible and non-extreme values. 
The parameterisation is probably conservatively small for 



very young stars that may be actively accreting (Hosokawa 
fc Omukai|[2009 i 



We introduced the intruder at a random orientation rel- 
ative to the hierarchy, with relative velocity at infinity v 
and the impact parameter distributed proportionally to the 
probability of the encounter (i.e. proportional to the square 
of the impact parameter) out to a maximum value. The max- 
imum value varied according to the stellar masses and the 
velocity of the encounter so that the maximum pericenter 
between the intruder and the hierarchy was equal to 3ai. 
If, at the end of an encounter, a stable triple still existed 
(according to the approximate triple stability criterion of 
Mardling fc Aarseth|20"oT I we sent another intruder n, ran- 



domly drawn in the same way as the first one. This process 
repeated until either the triple had been disrupted or else 
the system had survived 20 intruders. Throughout the paper, 
if a hierarchy was reduced to singles and binaries we refer 
to it as 'resolved'. 

For semi-major axes we used stable combinations of 
ao/au S {10,30,100} and oi/au G {100,300,1000}. With 
v/km s _1 G {1,5,20} we covered velocities appropriate to 
a wide range of cluster environments. For each set of hier- 
archy and encounter parameters we ran 16384 random re- 
alisations, for a total of nearly 3 x 10 5 individual systems, 
each undergoing up to 20 encounters. This limit was arbi- 
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Figure 2. The intruder and inner binary mass functions for all runs (thin lines) and runs where a collision occurred (thick lines) with ao 
= 10 au and ai = 300 au. Collisions occur more often with more massive intruders. The mass of the hierarchy components only becomes 
a factor when the hierarchy is no longer energetically hard with respect to the intruder. 



Table 1. Symbols describing the hierarchy and the intruder. 

mo Inner binary primary mass 

mi Inner binary secondary mass; mi <= mo 

m2 Outer binary secondary mass; mi <= mo 

Mb Inner binary mass; Mb = mo + mi 

Mt Triple mass; Mt = mo + mi + m2 

Mi Intruder mass 

ao Inner binary semi-major axis 

ai Outer binary semi-major axis 

v Relative velocity between the intruder and the triple 



trary; the number of encounters a hierarchy may experience 
depends sensitively on its environment. The simulations de- 
scribed in section [5] include the appropriate encounter rate 
for example clusters. 



3 SMALL-N RESULTS 

Table [3] summarises the collision frequency for each set of 
runs. For each setup we show the percent of runs with var- 
ious collision partners for all of the runs as well as for only 
those runs in which the original hierarchy resolved to a mix- 
ture of singles and lone binaries. We regard this quantity as 
more fundamental than the total percentage; the overall col- 
lision percentages can only rise as more and more intruders 
are simulated, while the resolved fraction is representative 
of the overall potential for collisions. The number of intrud- 
ers that might be expected in a real cluster depends on the 
timeframe of interest and the cluster properties. In section 
[5] we consider a limited set of these runs in a cluster setup 
similar to something like the Orion Nebula Cluster. We now 
discuss some of the trends apparent in these results; most of 
these are consistent with results from binary-single scatter- 
ing work. 



In figure [2] we show for one of the hierarchy setups 
(ao = 10 au, ai = 300 au; this is representative of all the 
sets) the intruder star Mi and inner binary Mb mass func- 
tions, showing both the input mass function and that for 
systems involved in collisions. At relatively low intruder ve- 
locities (v = 1 and 5 km s _1 ) the mass function of inner 
binaries in colliding systems is nearly identical to the input 
mass function. More massive intruders are more likely to 
result in collisions. At high relative velocities large values 
of Mt, likewise begin to become more likely. This is due to 
the increased energy of the intruders. Taking the hierarchy 
for the moment as a binary consisting of Mb and 7712, most 
systems are energetically hard to most intruders at 5 km 
s _1 with the large hierarchy masses we are dealing with. At 
20 km s~ a less massive hierarchy is not robust to ionisa- 
tion, and more massive systems are needed to remain intact 
for several collisions. 

When a collision does occur, the stars of the inner bi- 
nary are most likely to be involved. In figure [3] we show the 
fractional distribution of collision partners for the runs with 
ai = 10ao. Shades of blue are collisions involving at least 
one member of the inner binary, with the darkest shade in- 
dicating a collision between mo and mi . Collisions involving 
the outer member of the hierarchy are quite rare. 

In figure [4] we plot the collision fraction as a function 
of the outer to inner binary ratio, ai/ao- While we have 
too few points to determine the functional form of this re- 
lationship, the general trend is clear and unsurprising; the 
better the inner binary can be approximated as a single body 
the lower the collision fraction. The other relationship made 
clear in this plot is the decreasing collision fraction as the 
inner binary separation increases. This is analogous to the 
binary-single and binary-binary collision results of |Fregeau] 
et al. (20041, who show that the collision fraction increases 



roughly as a weak (sub-linear but not constant) power of 
R*/a. 
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Figure 3. For those runs with collisions, the distribution of colli- 
sion partners. Blue shades are double collisions involving a mem- 
ber of the inner binary. Red and green show the small number 
of collisions between the outer member of the hierarchy and the 
intruder and triple collisions, respectively. Collisions are predom- 
inantly between the members of the inner binary. 
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Figure 4. The fraction of runs that have a collision versus the 
ratio of the outer to inner binary semi-major axis. Thick lines 
show the result for all runs, thin lines show only resolved runs. 



Table 2. Collision fractions for the isolated encounter simu- 
lations. 

percent of runs with colliders 
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Semi-major axes ao and ai are in au; velocity v is in km 
s . Collision partners '0', T', and '2' are the primary, inner 
secondary, and outer secondary. Partner 'i' is the impactor. 
Each set shows the percentage of collisions for all runs as well 
as only resolved runs 



Figure [5] shows the collision fractions as a function of 
intruder velocity, again for the runs with ai — 10ao. The de- 
cline in the overall fraction is, as mentioned above, mainly 
due to the unresolved nature of an increasing percentage of 
runs as velocity increases. This tendency for higher velocity 
runs to require more encounters to resolve the hierarchy can 
be explained by appeal to the existing extensive work done 
on binary-single scattering ( Hut fc Bahcall|1983 Hut|1983 |. 
Treating the inner binary as its centre of mass particle (so 
that the hierarchy is reduced to a single binary with compo- 
nents Mb and ma), the critical energy at which the energy 
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Figure 5. The fraction of runs that have a collision for all the runs 
with ai = lOap as a function of the velocity of the intruder. Thick 
lines show the result for all runs, thin lines show only resolved 
runs. 



of the binary plus intruder is zero: 

v 2 = M b m 2 (M b + m 2 + Mj) 1 
'' ' Mi(M b + m 2 ) <zT 



(3) 



For roughly equal mass stars, this is within a factor of order 
unity of equation [2] For intruder velocities v < v c encoun- 
ters are likely to be more complex. As v gets higher relative 
to v c the most likely non-flyby outcome is ionisation of the 
binar y, but its cross secti on becomes increasingly small, as 
v~ 2 ( Hut fc Bahcall|l983| . As the timescale over which each 
individual encounter takes place shrinks, so does the chance 
of close interactions between the intruder and a binary com- 
ponent, and more encounters are flybys that leave the binary 
relatively unscathed. 

The decrease in the resolved fraction itself with increas- 
ing velocity makes sense in this approximation. The expla- 
nation for the increased collision fraction of resolved runs 
with increasing v has its roots in the same idea. In order for 
the hierarchy to be resolved by our definition, the possible 
outcomes are ionisation of the outer binary from the inner bi- 
nary, more complicated reorganisation of the four stars into 
binaries and singles, or a collision. As mentio ned above the 
ionisation cross section decreases as v~ 2 , while Fregeau et al. 



( 2004 1 showed that for binary-single and binary-binary en- 



counters the collision cross section drops to a constant but 
small value consistent with the physical cross section of the 
individual binary components. 

In this high v limit, the ionisation cross section can then 
drop below the collision cross section, particularly in a hi- 
erarchy when one member is a binary. The problem is then 
effectively a binary-single scattering problem with the bi- 
nary being the inner component of the hierarchy. That is, at 
high encounter velocities relative to the outer binary criti- 
cal velocity, the increased cross section of the hierarchy no 
longer helps to lower the inner binary's encounter timescale. 
Switching focus now to the inner binary and its critical ve- 
locity, with more massive stars it will have a lower value of 
v/vc, and thus higher collision rates ( jFregeau et al.|2004| > 

This is seen in figure [6] showing the fraction of colliding 
systems as a function of M b + Mi for the ai = 10ao runs 
for all three velocities (we plot these masses rather than the 
system total since the outer hierarchy member seldom took 



part in a collision). As v increases, the fraction of resolved 
runs drops, but those that are resolved are more likely to 
be resolved due to collisions than by disruption of the outer 
binary. The most massive systems have lower encounter ve- 
locities relative to the inner binary v c , and thus a higher 
collision fraction. For example, an inner 30 au binary with 
20 Mq components encountering a 5 Mq impactor at v — 20 
km s _1 has v/v c ~ 0.4. If we change the binary masses to 
2 Mq, and v/v c ~ 2.8. As v increases the problem reduces 
to an encounter with a single binary (the inner one), and 
the resolved collision fraction for different hierarchies start 
to converge, seen in figure [5] 

These arguments are clearly approximate, treating the 
hierarchy as disconnected binaries, but the results seem to 
fall in line with the extensive work done on cleaner binary- 
single and binary-binary encounters. The details of triple- 
single encounters have only recently begun to be explored 
in the same way ( |Leigh fc Geller|2012[ >, and more idealised 
experiments than those we have performed here (e.g. a more 
controlled mass function) are needed to tease out the funda- 
mental details. The more naturalistic setup we have chosen 
gives some hints about the trends that might be expected 
for the specific case of massive hierarchies in a young clus- 
ter. The question then becomes, how many hierarchies will 
be resolved in a realistic cluster, rather then the arbitrary 
20 encounters that we have performed? 



4 THE COPLANAR SPECIAL CASE 

The work we discussed above used hierarchies with ran- 
domly oriented orbital planes. While this is the most gen- 
eral case and appropriate to dynamically formed hierarchies, 
in a disc formation scenario nearly coplanar orbits might 
be expected. This special case caries with it two important 
physical effects. First, coplanar orbits maximise the net an- 
gular momentum of a hierarchy. Previous scattering work 
has shown that the angular momentum of a system is a key 
parameter in determining the scattering outcome (Mikkola 
fc Valtonen||1986 1. Specific to the situation at hand, Leigh 
& Geller (20121 showed that in triple-single encounters in- 



creasing the total angular momentum decreases the collision 
probability. 

Second, coplanar orbits suppress the Kozai-Lidov mech- 
anism ( |Kozai|[T962| |Lidov|[T962l |Naoz et al.pHl} in which 
secular perturbations can lead to extreme oscillations in the 
inner eccentricity of a inclined hierarchical systems like we 
study here. The timescale over which the oscillations take 
place is of order ( Kiseleva et al.|[l998 1 



TKL = 



_2_Ff_ 



(1 



mo + mi + m + 2 

I ) 

1712 



(4) 



with Po and eo the period and eccentricity of the inner bi- 
nary, and subscript 1 the outer binary. For the system masses 
and sizes we used, tkl can be of order 10 s yr. In general, the 
cumulative integration time of the small- N experiments was 
shorter than this, and the mechanism was already somewhat 
suppressed. Over the 2 Myr large- N runs that we discuss be- 
low, some of the systems would be able to go through one 
or more Kozai cycles, and potentially become tidally inter- 
acting or merge. 

As a topic of future work, the interplay between tidal 
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Figure 6. The fraction of runs with collisions versus the mass of the inner binary plus intruder for all the runs with a\ = 10ao- The 
shaded regions show the uncertainty in the estimate assuming Poisson counting statistics. The fractions are shown for all the runs done 
(black lines) and only resolved runs with no remaining triple systems (blue lines). Mean values over all bins are shown as the points at 
10 M Q . 



dissipation and the secular dynamics of hierarchical systems 
may be interesting in the case of massive primordial triples. 
The Kozai cycle with tidal friction (KCTF) mechanism has 
been extensively studied in lower-mass stars and planetary 
systems (e.g. |Mazeh & Shaham||1979| |Kiseleva et al.||1998 



Egglcton & Kisclcva-Egglcton 200 1[ |Fabrycky fc Tremaine] 



2007 



Perets & Naoz 2009 1. The short Kozai periods (relative 



to the main sequence lifetime of a massive star) and large 
radii of massive stars suggests that the process could be im- 
portant for producing massive short period binaries. This is 
not the point of our work here, however. In order to remove 
the possibility of merging inner binaries due solely to the 
internal dynamics of the hierarchy, we used co-planar hier- 
archies in the cluster experiments to remove the inclination- 
eccentricity exchange. 

In order to make a direct comparison between the clus- 
ter runs and the scattering results, we ran a limited set 
of scattering experiments directly tailored to the cluster 
characteristics that we describe fully below. While other- 



wise identical to the previous scattering results, these ones 
used coplanar hierarchies with v — 3. We performed 8096 
runs each of the following semi-major axis pairs: [oq,Oi] = 
[30, 300], [30, 1000], [10, 100], [10, 300] and [10,1000] au. Ta- 
ble [4] shows the results of these runs; the collision rate is 
about one half that of the general case, and the general 
trends with semi-major axes are similar. 



5 LARGE-N EXPERIMENTAL SETUP 

The small- N experiments we ran were necessarily highly ide- 
alised (a clustered environment was approximated by re- 
peatedly sending in single intruders). This was necessary in 
order to run a large number of experiments in a reasonable 
amount of time, and desirable in order to control the pa- 
rameters of the encounters. In order to verify the results in 
a more 'realistic' setting, we ran a limited set of full clus- 



ter experiments using the TV-body code NBODY6 (Aarseth 
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Table 3. Collision fractions for the coplanar isolated en- 
counter simulations. 

percent of runs with colliders 



«0 a \ v 
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total 


10 100 3 
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10 300 3 
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69% resolved 
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0.1 
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56% resolved 
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0.5 


1.5 
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77% resolved 
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12.6 
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2.6 
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0.2 
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59% resolved 
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0.2 
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0.3 


0.2 
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Semi-major axes ao and ai are in au; velocity v is in km 
s . Collision partners '0', '1', and '2' are the primary, inner 
secondary, and outer secondary. Partner 'i' is the impactor. 
Each set shows the percentage of collisions for all runs as well 
as only resolved runs 
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Figure 7. Symbols show the fraction of cluster runs that have a 
collision versus the outer semi- major axis ai (top), and the ratio 
of the hierarchy's outer to inner binary semi-major axis (bot- 
tom), plotted over the coplanar small- N scattering results shown 
as lines. The results at ai = 300 au have been offset horizontally 
for clarity. 



2000 



2003 1 . The code uses the 4th order Hermite algorithm 
( jMakino fc Aar seth 1992) to integrate the stars, augmented 
by algorithmic regularisation of close encounters and bina- 



ries ( Kustaanheimo & Stiefcl 196 5|. Stellar evolution is in- 
cluded via lookup tables ( jHurley et al.| 2000), providing the 
radii of the stars used in collision detection. 

We set up clusters of 2048 stars using the mass function 



described in |Maschberger (20121 between 0.1 and 1 50 Mq. 
This mass function is effectively the same as e.g. a |Kroupa| 
( 2001 1 or Chabrier ( 2003 1 mass function; we chose it mainly 



for computational ease. The mean mass of this mass func- 
tion is ~ 0.65 Mq, giving a cluster mass of ~ 1300 Mq, 



although the exact value was not fixed. We used King ( 1966 1 
models with Wo = 6 and a central mass density of 2 x 10 
Mq pc -3 . This translated to half- mass radii of ~ 0.4 pc, and 
one dimensional velocity dispersions of ~ 1.5 km s _1 . These 
characteristics were chosen to be broadly match the core 
properties of a young massive star forming region such as 
the Orion Nebula Cluster (Hillenbrand & Hartmann 19981. 
Since the encounter velocities should have been expected 
to be distributed roughly as a Maxwellian with dispersion 
v2~0"i w 2 km s~ and mode ~ 3 km s" 1 , these results should 
be directly comparable to the v — 3 coplanar small- N re- 
sults. 

We initialised the hierarchical systems in the follow- 
ing manner: after randomly selecting a star from those that 
were more massive than 10 Mq and inside roughly the half 
mass radius, we chose a binary mass ratio from the same 
power law as the small- N experiments, p(q) oc q~ 01 . We 
then adjusted the mass of the star that most closely matched 
this target mass to match it exactly, and moved it into 
orbit around the primary. We chose the third star's mass 
in the same way and placed it in a circular coplanar orbit 
around the inner binary's centre of mass. The hierarchy sep- 
arations we used were the same as the coplanar scattering 
results: [ao, ai] = [30, 300], [30, 1000], [10, 100], [10, 300] and 
[10, 1000] au. 

We ran 512 of these cluster simulations for each set 
of hierarchy separations, each for 2 Myr, and checked for 
collisions involving the hierarchy members. As controls we 
also ran 512 realisations each of identical setups but with 
only a 10, 30, 300 or 1000 au binary included. 



6 LARGE-N RESULTS 

As for the small- N experiments, for each run we checked for 
collisions involving members of the original hierarchy. If no 
collisions occurred we also checked if the hierarchical system 
was still intact by searching for hierarchies including at least 
one of the original inner binary members. This allowed us 
to determine both a collision fraction for the specific cluster 
setup that we simulated as well as an approximation to the 
'resolved only' results presented earlier. Table [4] shows the 
results for both the hierarchy runs and the control runs. In 
figure[7|we plot the collision fractions on top of the coplanar 
v = 3 km s _1 scattering results. The filled symbols show the 
fraction of collisions for all the cluster simulations, and the 
open symbols are the results for just the resolved runs. 

Across the full range of hierarchies we simulated, both 
the total and resolved collision fractions in the cluster runs 
are remarkably constant, with little evidence of the trends 
with semi-major axes seen in the scattering results. While 
there may be a hint of a lower collision fraction among re- 
solved runs as ai increases from 300 to 1000 au, the ao = 10, 
ai = 100 au runs had fewer collisions than their a± = 300 
au counterparts. We suspect the reason for this disagree- 
ment lies in the breakdown of the small- jV assumptions in 
the (simulated) reality of a cluster environment. With inter- 
action radii of order 1000 au, the mean interstellar separa- 
tion in the cores of the clusters were only a factor of a few 
larger than the encounter radius. The extremely short inter- 
action times there increased the number of stars that could 
be taking part in an interaction, and the encounters were no 
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longer well approximated as a four-body interaction. In the 
top panel of figure]?] note that the relative disagreement be- 
tween the scattering results and the cluster results increases 
with m, with only the smallest hierarchy (ai = 100 au) 
displaying agreement between the clusters and the isolated 
encounters. 

This situation is reminiscent of the work of Tanikawa 



10" 



et al. ( 2012 1, who showed that the actual formation of hard 
binaries in a core collapsed cluster involves many stars and 
bears little resemblance to the three-body approximation 
that is frequently used. Furthermore, with the very massive 
systems and low velocities we considered gravitational fo- 
cussing is dominant, and the impact parameters used in the 
scattering experiments were of the order of the size of the 
cluster core. It is not clear that the encounter geometries 
and characteristics in the cluster were well-described by the 
idealised scattering runs. 

For all of the setups, just over one-half of the hierarchies 
resolved down to binaries and singles over the 2 Myr inte- 
gration. If we had increased the integration time of the runs 
we expect that the collision fraction would have converge 
to the resolved-only result. In practice stellar evolution be- 
gins to complicate the interpretation of the runs past a few 
Myr, one of the reasons for our choice of a 2 Myr cutoff. The 
overall collision rates for all but the ao = 10, ai = 1000 au 
runs were roughly 2-3 times those of the binary-only runs. 
For the clusters (and integration times) that we simulated, 
the enhancement to collisions caused by the hierarchies was 
significant, but was a factor of order unity rather than an 
order of magnitude. 

Besides collisions, the heightened interaction that the 
inner binaries experience relative to their non-hierarchical 
counterparts should result in a higher frequency of stars es- 
caping the cluster. In figure [8] we plot the number of stars 
ejected per cluster with velocity greater than a given value, 
for escapers with mass greater than 1 Mq only. Our escaper 
selection was simple; stars at radii greater than 5 initial half 
mass radii and velocities v esc > 10 km s _1 were chosen (10 
km s _1 is about 5 times the escape velocity at 5 half mass 
radii, and safely outside the bulk of the cluster's velocity 
distribution; a star at that velocity would be unambiguously 
identified as a high velocity escaper). As expected the max- 
imum escape velocity increases with the binding energy of 
the binary for the binary-only runs, and with the binding 
energy of the inner binary for the hierarchies. Notably, the 
clusters hosting hierarchical systems eject high velocity es- 
capers at a rate 2-3 times greater than the binary-only 
clusters. 

While this result is limited (a single cluster configu- 
ration hosting a single hierarchy with two different archi- 
tectures), the issue may be worthy of more detailed study 
elsewhere. It suggests that detailed examinations of escaper 
quantities may be more sensitive to the initial conditions 
of massive stellar systems than has been appreciated. This 
issue is beginning to become more important as models of 
individual clusters grow more sophisticated, and variations 
in these quantities by factors of order unity are used to dis- 
criminate between different formation scenarios (e.g |Fujii| 
et al.|2012 i. 
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Figure 8. The mean number of stars per cluster with mass 
greater than 1 Mq that escape velocity v eS c greater than some 
value. Runs with initial binaries only are the lower curves in 
shades of gray. 



Table 4. Collision fractions for the full cluster simulations. 
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Semi-major axes ao and a\ are in au; velocity v is in km 
s -1 . Collision partners '0', T', and '2' are the primary, inner 
secondary, and outer secondary. Partner 'i' is the impactor. 
Each set shows the percentage of collisions for all runs as well 
as only resolved runs 

7 CONCLUSIONS 

After their formation wide, massive binaries, in typical 
young open clusters with modest velocity dispersions, will 
survive many encounters during their main sequence life- 
time. If the binary is instead a hierarchical triple, colli- 
sions may occur with interestingly high frequency over short 
timescales of only a few Myr. We carried out two series of 
numerical experiments to determine the collision rates of 
massive hierarchical star systems in young clusters. 

Our small- N scattering experiments suggest that if a 
hierarchical system is continuously bombarded by intruding 
stars until it resolves into binaries and singles, then tens of 
percent of the hierarchies will experience a collision. The 
original spacing of the hierarchy and the velocity of the in- 
truders determine the precise value. The question is then, 
how many hierarchies will get resolved into binaries and sin- 
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gles in a typical cluster environment? In our large- A simula- 
tions of a cluster similar to something like the Orion Nebula 
Cluster, roughly 50 per cent of our initial hierarchical sys- 
tems were broken up into binaries and singles in 2 Myr, with 
an overall collision fraction between about 7 and 14 per cent 
depending on the initial hierarchy configuration. For the spe- 
cific cluster we simulated, these numbers represent a mod- 
est increase to the collision fraction for lone binaries, which 
were under 5 per cent. The overall agreement between the 
scattering results and the full cluster simulations is at best 
suggestive. We suspect this is due to the breakdown of the 
scattering experiments' isolated-encounter assumption when 
the hierarchies are placed in a dense cluster; a lower-density 
environment simulated over longer time periods might pro- 
vide a better match, but then the issue of the massive stars' 
short lifetimes rears its head. 

In addition to collisions, the high rate of encounters 
produces a similar excess of runaway stars compared to sys- 
tems with no primordial hierarchies. The relative scarcity of 
collision products and high velocity escapers produced by 
a cluster can make them useful as a diagnostic of different 
cluster initial conditions or formation scenarios. Our results 
suggest that the higher order multiplicity characteristics of 
massive stars are an important input to any dynamical ex- 
periments in which they play a role in producing these rare 
events. 
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